Integration

What if we wanted to look at population in 20 years given an initial condition

Two options

• explicit solution to differential equation is known; e.g. you can integrate both sides of the equation! Not always possible but lets look at a case where it is possible

• must be solved by iteration; this is what we do when we can’t integrate both sides

Explicit Solution is available

source("../R/exppop.R")

exppop

## function (T, P0, r, K) 
## {
##     P = P0 * exp(r * T)
##     if (P > K) {
##         P = K
##     }
##     return(P)
## }

# gives population after any time given an initial population

# 20 rabbits, growth rate of 0.01 how many in 30 years
exppop(T=30, P0=20, r=0.01, K=1000)

## [1] 26.99718

# if we want to see how population evolves over time - generate a time series by running our model for each point in time

initialrabbits = 20
years = seq(from=1, to=100, by=2)
Ptime = years %>% map_dbl(~exppop( P0=initialrabbits, r=0.01, K=1000, T=.x))

# keep track of what times we ran
Ptime = data.frame(P=Ptime, years=years)

ggplot(Ptime, aes(years,P))+geom_point()+labs(x="years",y="Rabbit Population")

Difference Equations

Population models can be discrete (rather than continuous)

So we could implement them as difference equations and iterate

source("../R/discrete_logistic_pop.R")
# notice how a for loop is used to iterate

# how many rabbits after 50 years given a growth of 0.1
# starting with 1 rabbit - but a carrying capcity of 500

discrete_logistic_pop

## function (P0, r, K, T = 10) 
## {
##     pop = P0
##     for (i in 1:T) {
##         pop = pop + r * pop
##         pop = ifelse(pop > K, K, pop)
##     }
##     return(pop)
## }

discrete_logistic_pop(P0=1, r=0.05, K=200, T=50)

## [1] 11.4674

# save results
discrete_result = discrete_logistic_pop(P0=1, r=0.05, K=200, T=50)

# lets also keep the parameters for use later
P0=1
r=0.05
K=200
T=50

Compare discrete and analytic results

source("../R/exppop.R")

exppop(P0=P0, r=r, K=K, T=T)

## [1] 12.18249

analytic_result = exppop(P0=P0, r=r, K=K, T=T)

analytic_result

## [1] 12.18249

discrete_result

## [1] 11.4674

# why are they different
# look at trajectories

growth_result = data.frame(time=seq(from=1,to=100))

growth_result$Panalytic = growth_result$time %>% map_dbl(~exppop( P0=1,r=0.05, K=200,T=.x ))

growth_result$Pdiscrete = growth_result$time %>% map_dbl(~discrete_logistic_pop( P0=1,r=0.05, K=200,T=.x ))

tmp = growth_result %>% gather(key="Ptype",value="P",-time)
ggplot(tmp, aes(time,P, col=Ptype))+geom_point()

# try running them for longer time periods to see what happens 
# change the value of r, K , P0 - see how it effects the results

Solving using numeric integration

Using a solver….when you can’t do the integration by hand :)

For example, if you added a carrying capacity as threshold where it stops growing

In this case we integrate by iteration - approximates analytic integration

Numerical integration with ODE

Implement the differential equation as a function that

• returns the derivative (as a list)

*inputs time, the variable(s) and a parameter list

(it needs time even though you don’t use it)

My convention: name derivative functions starting with d to remind myself that they are computing a derivative

ODE

Only works for Ordinary Differential Equations - single independent variable (in our case time)

Partial differential equations - more than 1 independent variable (e.g x and y if changing in space)

R has a solver called ODE for solving ordinary differential equations from package desolve

ODE requires

• initial conditions
• values of independent where we want values of dependent variable (e.g times where we want population)
• the derivative as a function
• a list that contains all parameter values (or if you have only one parameter then you can use a single value)

source("../R/dexppop.R")

dexppop

## function (time, P, r) 
## {
##     dexpop = r * P
##     return(list(dexpop))
## }

library(deSolve)
initialrabbits = 20
years = seq(from=1, to=100, by=2)

# run the solver
Ptime = ode(y=initialrabbits, times=years, func=dexppop,parms=c(0.01))
head(Ptime)

##      time        1
## [1,]    1 20.00000
## [2,]    3 20.40404
## [3,]    5 20.81623
## [4,]    7 21.23675
## [5,]    9 21.66576
## [6,]   11 22.10344

colnames(Ptime)=c("year","P")

# notice that there are additional pieces of information year, including the method used for integration
attributes(Ptime)

## $dim
## [1] 50  2
## 
## $dimnames
## $dimnames[[1]]
## NULL
## 
## $dimnames[[2]]
## [1] "year" "P"   
## 
## 
## $istate
##  [1]   2  52 105  NA   6   6   0  36  21  NA  NA  NA  NA   0   1   1  NA  NA  NA
## [20]  NA  NA
## 
## $rstate
## [1]  2.00000  2.00000 99.98839  0.00000  1.00000
## 
## $lengthvar
## [1] 1
## 
## $class
## [1] "deSolve" "matrix" 
## 
## $type
## [1] "lsoda"

# this also means you need to extract just the data frame for plotting
ggplot(as.data.frame(Ptime), aes(year,P))+geom_point()+labs(y="Population", "years")

# this also works (of course function can be by order)
Ptime=ode(initialrabbits, years, dexppop,0.03)
colnames(Ptime)=c("year","P")
ggplot(as.data.frame(Ptime), aes(year,P))+geom_point()+labs(y="Population", "years")

additional parameters

You can play a bit with changing your function to something that you can’t integrate “by hand”

BUT we might want more parameters

to work with ODE, parameters must all be input as a single list; similar to how we return multiple outputs from a function

see example below..lets add a carrying capacity

R code with carrying capacity

source("../R/dexppop_play.R")

dexppop_play

## function (time, P, parms) 
## {
##     dexpop = parms$r * P
##     dexpop = ifelse(P > parms$carry_capacity, 0, dexpop)
##     return(list(dexpop))
## }

# create parameter list
initalrabbits=2
newparms = list(r=0.03, carry_capacity=300)

#apply solver
results=ode(initialrabbits, years, dexppop_play,newparms)
head(results)

##      time        1
## [1,]    1 20.00000
## [2,]    3 21.23675
## [3,]    5 22.54993
## [4,]    7 23.94435
## [5,]    9 25.42499
## [6,]   11 26.99719

# add more meaningful names
colnames(results)=c("year","P")

# plot
ggplot(as.data.frame(results), aes(year,P))+geom_point()+labs(y="Population", "years")

# try again with different parameters
alternativeparms = list(r=0.04, carry_capacity=500)
results2=ode(initialrabbits, years, dexppop_play,alternativeparms)

## DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
##       taken on this call before reaching TOUT     
## In above message, I1 = 5000
##  
## In above message, R1 = 81.5108
## 

## Warning in lsoda(y, times, func, parms, ...): an excessive amount of work (>
## maxsteps ) was done, but integration was not successful - increase maxsteps

## Warning in lsoda(y, times, func, parms, ...): Returning early. Results are
## accurate, as far as they go

# look at results
head(results2)

##      time        1
## [1,]    1 20.00000
## [2,]    3 21.66575
## [3,]    5 23.47022
## [4,]    7 25.42499
## [5,]    9 27.54257
## [6,]   11 29.83650

colnames(results2)=c("year","P_parm2")

# plot
ggplot(as.data.frame(results2), aes(year,P_parm2))+geom_point()+labs(y="Population", "years")

# compare by combining into a single data frame
both = inner_join(as.data.frame(results), as.data.frame(results2))

## Joining, by = "year"

both_p = both %>% gather(key="model", value="P", -year)
ggplot(both_p, aes(year,P, col=model))+geom_point()+labs(y="Population", "years")

# try playing on your own - modify the function in some way

Differential Equation, Difference (Iteration by hand) comparison

Remember we have 3 ways now to calculate population

analytical solution - based on integration (exppop.R) BEST

using an ode solver for numerical approximation (exppop_play.R)

numerical integration using in discrete steps (discrete_logistic_pop.R)

Compare analytical, difference and ODE

Finally look at continuous derivative using ODE solve Needs initial conditions differential equation *parameters

source("../R/dexppop_play.R")

dexppop_play

## function (time, P, parms) 
## {
##     dexpop = parms$r * P
##     dexpop = ifelse(P > parms$carry_capacity, 0, dexpop)
##     return(list(dexpop))
## }

# set up using the same parameters
pcompare = list(r=r,carry_capacity=K)


# now run our ODE solver
result = ode(P0, growth_result$time, dexppop_play, pcompare)
head(result)

##      time        1
## [1,]    1 1.000000
## [2,]    2 1.051273
## [3,]    3 1.105172
## [4,]    4 1.161835
## [5,]    5 1.221404
## [6,]    6 1.284027

# we already have time - so just extract population
growth_result$Pdifferential = result[,2]

# compare all 3 approaches
tmp = growth_result %>% pivot_longer(cols=-time,names_to="Ptype",values_to="P")
ggplot(tmp, aes(time,P, col=Ptype))+geom_point()

# notice Pdifferential is closer to Panalytic than Pdiscrete

Other Examples

All differential and difference equations are approximations The analytical solution is exact

Notice that differential equations is a bit more accurate!

Diffusion

Diffusion can be implemented as a partial differential equation

Complicated to solve - but there are tool in R for specific types of partial differential equations

More info on differential equations in R

Diffusionn would require partial derivatives - time and space! it gets much more tricky …beyond this course

But we can appoximate diffusion with a difference equation - and iterative to get an estimate of how diffusion works

Example of Diffusion - difference equation implementation to see what some issues can be

Diffusion in Theory

Diffusion in 1 dimension

Diffusion in one dimension through time

Diffusion data structure

R implementation

source("../R/diffusion.R")

# run our diffusion model (iterative difference equation) with initial concentration of 10, for 8 timestep (size 1m), and 10 space steps (size 1s)
# using diffusion parameters 0.5 s/m2, 10 m2
result = diff1(initialC=10, nx=10, dx=1, nt=8, dt=1, D=0.5, area=10)

# a list is returned with our 3 data frames for concentration (conc), qin and qout
result

## $conc
##           [,1]     [,2]     [,3]      [,4]      [,5]        [,6]        [,7]
## [1,] 10.000000 0.000000 0.000000 0.0000000 0.0000000 0.000000000 0.000000000
## [2,]  7.500000 2.500000 0.000000 0.0000000 0.0000000 0.000000000 0.000000000
## [3,]  6.250000 3.125000 0.625000 0.0000000 0.0000000 0.000000000 0.000000000
## [4,]  5.468750 3.281250 1.093750 0.1562500 0.0000000 0.000000000 0.000000000
## [5,]  4.921875 3.281250 1.406250 0.3515625 0.0390625 0.000000000 0.000000000
## [6,]  4.511719 3.222656 1.611328 0.5371094 0.1074219 0.009765625 0.000000000
## [7,]  4.189453 3.142090 1.745605 0.6982422 0.1904297 0.031738281 0.002441406
## [8,]  3.927612 3.054810 1.832886 0.8331299 0.2777100 0.064086914 0.009155273
##              [,8] [,9] [,10]
## [1,] 0.0000000000    0     0
## [2,] 0.0000000000    0     0
## [3,] 0.0000000000    0     0
## [4,] 0.0000000000    0     0
## [5,] 0.0000000000    0     0
## [6,] 0.0000000000    0     0
## [7,] 0.0000000000    0     0
## [8,] 0.0006103516    0     0
## 
## $qout
##           [,1]     [,2]     [,3]     [,4]       [,5]       [,6]        [,7]
## [1,] 25.000000 0.000000 0.000000 0.000000 0.00000000 0.00000000 0.000000000
## [2,] 12.500000 6.250000 0.000000 0.000000 0.00000000 0.00000000 0.000000000
## [3,]  7.812500 6.250000 1.562500 0.000000 0.00000000 0.00000000 0.000000000
## [4,]  5.468750 5.468750 2.343750 0.390625 0.00000000 0.00000000 0.000000000
## [5,]  4.101562 4.687500 2.636719 0.781250 0.09765625 0.00000000 0.000000000
## [6,]  3.222656 4.028320 2.685547 1.074219 0.24414062 0.02441406 0.000000000
## [7,]  2.618408 3.491211 2.618408 1.269531 0.39672852 0.07324219 0.006103516
## [8,]  0.000000 0.000000 0.000000 0.000000 0.00000000 0.00000000 0.000000000
##      [,8] [,9] [,10]
## [1,]    0    0     0
## [2,]    0    0     0
## [3,]    0    0     0
## [4,]    0    0     0
## [5,]    0    0     0
## [6,]    0    0     0
## [7,]    0    0     0
## [8,]    0    0     0
## 
## $qin
##      [,1]      [,2]     [,3]     [,4]     [,5]       [,6]       [,7]
## [1,]    0 25.000000 0.000000 0.000000 0.000000 0.00000000 0.00000000
## [2,]    0 12.500000 6.250000 0.000000 0.000000 0.00000000 0.00000000
## [3,]    0  7.812500 6.250000 1.562500 0.000000 0.00000000 0.00000000
## [4,]    0  5.468750 5.468750 2.343750 0.390625 0.00000000 0.00000000
## [5,]    0  4.101562 4.687500 2.636719 0.781250 0.09765625 0.00000000
## [6,]    0  3.222656 4.028320 2.685547 1.074219 0.24414062 0.02441406
## [7,]    0  2.618408 3.491211 2.618408 1.269531 0.39672852 0.07324219
## [8,]    0  0.000000 0.000000 0.000000 0.000000 0.00000000 0.00000000
##             [,8] [,9] [,10]
## [1,] 0.000000000    0     0
## [2,] 0.000000000    0     0
## [3,] 0.000000000    0     0
## [4,] 0.000000000    0     0
## [5,] 0.000000000    0     0
## [6,] 0.000000000    0     0
## [7,] 0.006103516    0     0
## [8,] 0.000000000    0     0

# used filled contour to plot results
head(result$conc)

##           [,1]     [,2]     [,3]      [,4]      [,5]        [,6] [,7] [,8] [,9]
## [1,] 10.000000 0.000000 0.000000 0.0000000 0.0000000 0.000000000    0    0    0
## [2,]  7.500000 2.500000 0.000000 0.0000000 0.0000000 0.000000000    0    0    0
## [3,]  6.250000 3.125000 0.625000 0.0000000 0.0000000 0.000000000    0    0    0
## [4,]  5.468750 3.281250 1.093750 0.1562500 0.0000000 0.000000000    0    0    0
## [5,]  4.921875 3.281250 1.406250 0.3515625 0.0390625 0.000000000    0    0    0
## [6,]  4.511719 3.222656 1.611328 0.5371094 0.1074219 0.009765625    0    0    0
##      [,10]
## [1,]     0
## [2,]     0
## [3,]     0
## [4,]     0
## [5,]     0
## [6,]     0

filled.contour(result$conc, xlab="Time", ylab="Distance")

# or if you prefer this orientation (Distance on x axis)
filled.contour(t(result$conc), ylab="Time", xlab="Distance")

Change parameters (diffusivity D, and space and time steps (dx, dt))

# changes diffusivity and other parameters particularly
# diffusivity, dx and dt
res=diff1(initialC=10,nx=10,dx=1,nt=10,dt=30,D=0.006,area=1)

filled.contour(res$conc, xlab="Time", ylab="Distance")

# we can also see how much material moved from place to place each time step
filled.contour(res$qin, xlab="Time", ylab="Distance")

# play with time step, space step and parameters

Play

Try running the diffusion model with different time steps, space steps and parameters

Some examples with different parameters and space/time steps

# what if we increase diffusivity
resfast=diff1(initialC=10,nx=10,dx=0.5,nt=10,dt=10,D=0.08,area=1)
filled.contour(resfast$conc, xlab="Time", ylab="Distance")

# Discretization Issue Example
resunstable=diff1(initialC=10,nx=10,dx=1,nt=10,dt=10,D=0.8,area=1)
filled.contour(resunstable$conc, xlab="Time (fraction of hour)",ylab="Distance Along Path (m)", main="Pollutant Diffusion")

# this illustrates the problem with difference equations (and the challenges that methods for numerical integration try to overcome)
# if things are changing quickly we need to use much smaller time, space steps to avoid overshoot and instability

# so lets cut our step size by 10 (dt) (but then add 10 more steps (nx) to cover the same distance)
resunstable=diff1(initialC=10,nx=100,dx=1,nt=10,dt=1,D=0.8,area=1)
filled.contour(resunstable$conc, xlab="time",ylab="Distance Along Path")